Perturbation detection

ABSTRACT

A method for detecting perturbation of a physical system from a reference state associated with a reference parameter (ω 0 ) to a perturbed state associated with a perturbed parameter (ω) includes firstly deriving the reference parameter (ω 0 ). A reference vector (F) is then derived which describes the system&#39;s state at the reference parameter (ω 0 ). A measurement-related vector (Z) associated with a perturbed state of the system is then subtracted from the reference vector (F) to provide an error vector (E). The error vector members are summed and normalised by division of a summation of elements of a vector (F′) representing a derivative (f(ω 0 δω) e ) of a reference itself represented by the reference vector (F), the derivative (f(ω 0 δω) e ) being evaluated at the reference parameter (ω 0 ). This provides a result equal to the difference (ω−(ω 0 ) between the perturbed parameter (ω) and the reference parameter (ω 0 ).

BACKGROUND OF THE INVENTION

(1) Field of the Invention

This invention relates to detecting perturbation, and more particularlyto a method, an apparatus and a computer programme for detectingperturbation of a physical system's state from a reference state. Auseful application of this invention is to detection of perturbation ofthe state of a sensitive detection system such as a porous siliconbiosensor by an optical technique.

(2) Description of the Art

Porous silicon is known to be suitable for fabrication into sensitivedetectors. The nature of porous silicon allows a ready response at amicroscopic level to perturbations caused by environmental changes.Moreover, porous silicon lends itself to interferometric techniques:well-resolved Fabry-Perot fringes may be generated on illumination of alayer of the material providing a basis for sensitive interferometricdetection of its microscopic response. The response may simply beabsorption or trapping of molecules in the pores of the silicon or itmay be made more specific by grafting a recognition agent to an internalsurface of the pores. In either case, any molecular event occurringwithin a porous silicon layer will affect the layer's opticalproperties, resulting in a change in an interference fringe patternwhich can be detected and quantified if of sufficient magnitude.

Sailor et al. in WO99/12305 describe a porous silicon-based biosensor.Biosensors in general consist of two components: a recognition agentwhich reacts to produce a molecular or chemical signal in the presenceof a specific reagent and a signal transducer which converts themolecular recognition event into a quantifiable signal. The recognitionagent used in a porous silicon biosensor may be, for example, anantibody grafted into the layer and so a molecular reaction will occurin the presence of a specific antigen. Occurrence of such a reaction isobserved via an optical interference pattern.

In any periodic interference pattern, such as produced by the biosensordescribed in WO99/12305, the separation between peaks is proportional tothe optical path difference between alternative routes taken to adetector as the incident beam is partially reflected at interfaces.Under such circumstances, an accepted analysis technique is to locatepeaks in the Fourier transform spectrum of the interference pattern. Thelocation of such peak(s) provides an indication of fringe spacing(s) inthe original pattern, from which a value for optical path difference canbe derived. As will be described in more detail later, the interferencepattern developed using a porous silicon biosensor is collected as adiscrete set of spectral data points. One of the standard fast Fouriertransform (FFT) algorithms is therefore used as the basis for Fourieranalysis. The details of such an approach to analysing signals from aporous silicon biosensor are described by Janshoff A. et al. in J. Am.Chem. Soc. 120 pp 12 108-12 116 (1998). There are however variousdrawbacks to using this method of analysis, particularly in situationssuch as provided by the porous silicon biosensor in which the monitoredbinding event is capable of highly sensitive discrimination. Althoughthe Fourier transform method per se (i.e. finding the position of thepeak in amplitude of the Fourier transform of the data) is known to beoptimal in locating the frequency of a single tone in white noise, manyactually measured signals do not sufficiently approach this ideal and sodo not allow full exploitation of the method's capabilities. Moreover,use of the FFT in particular involves interpolation in the transformspectrum between regular sampling points. This severely limits theaccuracy with which a peak in the transform can be located.

There is generally a need for increased sensitivity in any detectionsystem, and this can be achieved either by rendering a detecting mediummore responsive to events to be detected or by improving extraction ofrequired data from an inevitably-noisy experimentally detected signal.More sensitive biosensors will aid, for example, early detection of adisease antigen or rapid detection of the presence of explosives. It isan object of the present invention to provide an alternative form ofperturbation detection.

SUMMARY OF THE INVENTION

The present invention provides a method of detecting perturbation of aphysical system from a reference state associated with a referenceparameter to a perturbed state associated with a perturbed parameter(ω), characterised in that the method includes the steps of:

-   -   a) deriving the reference parameter;    -   b) deriving a reference entity which describes the system's        state at the reference parameter;    -   c) deriving an error entity from the difference between the        reference entity and a measurement-related entity associated        with a perturbed state of the system; and    -   d) using the error entity to provide an indication of whether or        not the system state has been perturbed from the state        associated with the reference parameter, with the perturbed        parameter having become unequal to the reference parameter.

The expression ‘entity’ is used to mean a vector member, a vector or afunction.

The invention provides for detection of perturbation with respect to areference parameter and responds to a change from this indicated by theperturbation parameter. Such change is normally small, and so theinvention provides the advantage that it makes it possible to employmathematical approximations which are applicable to small changes inparameters.

The error entity may be a vector with multiple members and step d) thenincludes calculating the sum of the error vector's members to obtain anindication of whether or not the perturbed parameter has become unequalto the reference parameter. The error entity may have memberscharacterised by relatively high signal to noise ratio compared otherpossible members which might otherwise be selected for deriving it. Itmay have members derived from a region of a peak in a derivative of thereference entity.

Prior to derivation of the error entity in step c), themeasurement-related entity is normalised by projection on to a spaceorthogonal to that of the reference entity. Step d) may includedetermining the difference between the perturbed and referenceparameters by error entity normalisation with respect to an entity whichis the summation of elements of an entity representing a derivative of areference function evaluated at the reference parameter, the referencefunction being represented by the reference entity.

Step b) may incorporate:

-   -   a) predicting from the reference parameter a position of a peak        in a Fourier transform spectrum of observation data;    -   b) selecting a subset of the observation data over the predicted        peak calculating a direct Fourier transform (as herein defined)        for the subset; and    -   c) analysing the direct Fourier transform in order to derive a        more accurate determination of the position of the peak.

The reference entity may be derived in step b) by a process whichincludes filtering by projection of an entity on to a set ofpre-determined entities having a range of arguments all of which differfrom the reference parameter by less than one tenth of the referenceparameter.

Step c) of the method may comprise calculating the error entity by aprocess including deriving the measurement-related entity by a processwhich includes data filtering by projection of a measurement entity onto the set of pre-determined entities.

The reference entity may be derived in step b) from an average of aseries of observation data vectors. These vectors may be time-varyingand this may be taken into account when their average is derived. Thereference entity may be a sum of a static reference entity and alinearly time-varying entity.

The physical system may be a sensor pervaded by a medium having variablecomposition. The sensor may be a porous silicon sensor with porespervaded by a solvent medium, the perturbed and reference parameters arethen optical thicknesses of a region of the sensor pervaded by differentsolvent compositions and in which interference patterns are generatedfor derivation of the reference entity and error entity.

In an alternative aspect, the present invention provides apparatus fordetecting perturbation of a physical system from a reference stateassociated with a reference parameter to a perturbed state associatedwith a perturbed parameter, characterised in that the apparatusincludes:

-   -   a) a monitoring device for monitoring the physical system to        obtain indications of the system's state;    -   b) a computer system programmed to process the system state        indications to:        -   i) derive the reference parameter (ω₀);        -   ii) derive a reference entity (F) which describes the            system's state at the reference parameter (ω₀);        -   iii) derive an error entity (E) from the difference between            the reference entity (F) and a measurement-related            entity (Z) associated with a perturbed state of the system;            and        -   iv) use the error entity (E) to provide an indication of            whether or not the system state has been perturbed from the            state associated with the reference parameter (ω₀), with the            perturbed parameter (ω) having become unequal to the            reference parameter (ω₀).

In a further alternative aspect, the present invention provides computersoftware for use in detecting perturbation of a physical system from areference state associated with a reference parameter to a perturbedstate associated with a perturbed parameter, characterised in that thesoftware contains instructions for controlling a computer to implementthe steps of:

-   -   a) deriving the reference parameter;    -   b) deriving a reference entity which describes the system's        state at the reference parameter;    -   c) deriving an error entity from the difference between the        reference entity and a measurement-related entity associated        with a perturbed state of the system; and    -   d) using the error entity to provide an indication of whether or        not the system state has been perturbed from the state        associated with the reference parameter, with the perturbed        parameter having become unequal to the reference parameter.

The apparatus and software aspects of the invention may have preferredfeatures equivalent mutatis mutandis to those of the method aspect.

DESCRIPTION OF THE FIGURES

Embodiments of the invention will now be described, by way of exampleonly, with reference to the accompanying drawings, in which:

FIG. 1 is a schematic illustration of a porous silicon biosensor for usein perturbation detection in accordance with the invention;

FIG. 2 is a typical interference pattern obtained with the biosensor ofFIG. 1;

FIG. 3 is a spectrum obtained from the interference pattern of FIG. 2using a prior art FFT analysis technique;

FIG. 4 is a schematic illustration of an idealised result of monitoringoptical path difference over time, this result being one whichtheoretically may be obtained after analysis of a changing interferencepattern obtained using the biosensor of FIG. 1;

FIG. 5 shows three plots equivalent to that shown in FIG. 4, butextracted by a prior art signal analysis technique from an actual dataset obtained using the biosensor illustrated in FIG. 1;

FIG. 6 is a flow chart of steps involved in deriving a filteredcalibration function for use in the present invention;

FIG. 7 illustrates spectra equivalent to those shown in FIG. 5, thistime obtained using a first stage of the present invention;

FIG. 8 is a rescaled version of the plots of FIG. 7;

FIG. 9 is a flow chart of the steps involved in analysing data inaccordance with the analysis system of the present invention;

FIG. 10 provides graphs of an estimated reference function ƒ(ω₀, k)_(e)and a cos(ωk) approximation to it plotted against wavenumber in eachcase.

FIG. 11 provides graphs of error functions plotted against wavelength.

FIG. 12 shows regions of the FIG. 11 graphs at a peak of one of thegraphs;

FIG. 13 illustrates spectra equivalent to those plotted in FIGS. 5 and7, except that two stages of processing in accordance with the inventionwere employed;

FIG. 14 shows measured and calculated interference patterns obtainedusing porous silicon sensor with eighteen layers and forming a Braggmirror; and

FIG. 15 is a graph of response of the multilayer sensor processingsystem used to obtain the patterns shown in FIG. 14, and is plottedagainst time sample number.

DESCRIPTION OF A PREFERRED EMBODIMENT

Referring to FIG. 1 a porous silicon biosensor 10 can be used to provideexperimental results which are particularly suited for perturbationdetection by embodiments of the present invention. Basic principles ofuse of the biosensor 10 are as follows. The biosensor 10 comprises aporous silicon layer 12 on a bulk (substantially non-porous) siliconsubstrate 14. When the porous layer 12 is illuminated by light 16 from awhite light source 18, the light 16 is reflected at interfaces 20, 22between the porous silicon layer 12 and air and between the poroussilicon layer 12 and the bulk substrate 14. Reflected light is detectedby a CCD detector 24. Reflection at the two interfaces 20, 22 gives riseto two reflected beams which interact with each other to generate aninterference pattern 26, which is detected at the detector 24. Signalsdetected by the detector 24 are processed by a computer system 27.

The interference pattern 26 is essentially a measurement of theintensity of light reflected from the biosensor 10 as the lightwavelength is varied across the white light spectrum. In practice,alternative measures of intensity such as reflectance (in the pattern 26shown in FIG. 1) or voltage developed at the detector 24 (see FIG. 2)are used to provide an indication of the amplitude of the interferencepattern 26. Whatever measure is used however, interference peaks willoccur at wavelengths λ_(max) when the well-known interference conditionis satisfied:2n(λ_(max))d=mλ _(max)  (1)where d is the thickness of the porous silicon layer 12, n(λ_(max)) isthe refractive index of the layer 12 at wavelength λ_(max), is aninteger and normal incidence of the light 16 on the layer 12 is assumed.

A magnified view 28 of a region of the porous silicon layer 12 is alsoshown in FIG. 1. This view 28 depicts in microstructure a representationof the porous silicon material 30 on to which is grafted a recognitionagent 32. In this instance, the recognition agent 32 is illustratedbound to an analyte 34, for which it is a specific identifier. Therecognition agent 32 might, for example, be an antibody for a chemicalagent or for an explosive.

In the absence of the analyte 34, the porous silicon 12 with graftedrecognition agent 32 exhibits particular optical properties, and theseproperties enable an interference pattern to be generated with peaks inaccordance with Equation (1). If the analyte 34 is present in theenvironment however, then some proportion will bind to the recognitionagents 32 within the porous layer 12 which in turn will affect theoptical properties of that layer. Specifically, layer thickness willremain unchanged, but the refractive index n(λ_(max)) of the layer willbe affected by the additional bonded molecule. This will be evidenced asa movement in the positions of the interference peaks, as dictated byEquation (1).

It is to be noted that, strictly, the refractive index n of the layer iswavelength dependent, hence use of the n(λ_(max)) notation. In apractical application of the method described here, it is assumed thatthe layer refractive index is independent of wavelength and use is madeof an average value over the range of illuminating wavelengths employed.Accordingly, the layer refractive index referred to in Equation (1) willbe denoted n_(avg) when descriptions of practical implementations aregiven.

It is readily apparent from Equation (1) that the spacing betweenadjacent maxima in the interference pattern provides an indication ofn_(avg)d. If the interference spectrum in the absence of the analyte isknown, then observation of the interference spectrum 26 developed from abiosensor 10 in an unknown environment will provide an indication as towhether the analyte 34 is present or absent in that environment. Asinterferometry may be extremely sensitive, it is to be expected that, inprinciple, even small amounts of bound analyte can be detected usingthis biosensor. Firstly however, some means of analysing andinterpreting the interference spectrum must be developed.

FIG. 2 illustrates a typical example of an interference pattern 50generated from the biosensor 10 of FIG. 1. At various wavelengths oflight, the voltage developed at the detector 24 is measured. Suchmeasurements contribute to an interference pattern 50 comprising avector of voltage samples taken at a series of selected wavelengths.Janshoff et al. (referred to earlier) describe one approach toextracting a value of n_(avg)d from such an interference pattern 50. Themeasured data is first converted to represent a plot of the interferencespectrum against the reciprocal of the wavelength (equivalently,wavenumber k) and then the Fourier transform is taken. This transformsthe data to a “frequency” domain (in this context use of the word“frequency” refers to the frequency of the peaks in the spectralinterference pattern, and not to the frequency of the light used togenerate it). A peak in amplitude of the Fourier transform gives a valuefor the period of the peaks in the original pattern. Janshoff et al.describe use of a fast Fourier transform (FFT) method to generate aninterference frequency spectrum from the acquired data points. Followingthis method, an FFT spectrum 60 of the data shown in FIG. 2 isillustrated in FIG. 3. Two peaks 62, 64 are visible in the transformeddata. The larger peak 62 in the region of 0 nm separation values is dueto the zero frequency component of the interferogram 50 i.e. the offsetapparent in FIG. 2. The smaller peak 64 in the region of 9,500 nmcorresponds to a value of n_(avg)d, where n_(avg) is thewavelength-averaged value of n(λ_(max)), which gives rise to thespectral interferogram 50. In addition there is high frequency noise,visible as oscillations superimposed on the interferogram 50 of FIG. 2,but their contribution to the Fourier transform is not visible on thescale of FIG. 3. Sampled data points 66 a, 66 b, 66 c, etc., are alsoindicated in the plot of FIG. 3. These represent regular sampling pointsat which the value of the FFT function is calculated. In order to locatethe peaks, the interferogram data must be interpolated from the set ofwavenumbers k at which detector measurements are made.

Using this approach n_(avg)d can be measured to about 0.5 nm. For atypical layer thickness of 3000 nm, n_(avg)d is about 4000 nm, and sothis represents an accuracy of about 1 part in 10 000.

FIG. 4 illustrates schematically an idealised representation of a set ofresults which may be obtained using a chemical sensor variant of thebiosensor 10 of FIG. 1 and theoretically efficient signal processing.The results are displayed as a plot 70 of the value of n_(avg)d againstthe time at which each respective spectral interferogram 50 wascollected. These idealised readings have been assumed to have been takenwith a chemical sensor which does not contain the recognition agent 32of the biosensor 10. Instead certain molecules introduced into theenvironment will be simply trapped in the pores of the porous siliconlayer. In this, and experiments obtaining actual data which will bereferred to later, water is pumped through the sensor 10, which presentsa thin film to the flow. Initially the water is pure water, then asucrose solution is introduced before a return to the flow of purewater. When sucrose molecules are present, they are absorbed into theporous silicon layer.

The plot 70 comprises a series of data points 72 which, on average,follow a sloping background 74. The data points 72 show a clear rise inthe value of n_(avg)d at time t₁ and a clear drop in its value to thelevel of the sloping background 74 at t₂. The rise at t₁ corresponds tothe time at which sucrose was introduced into the water being pumpedthrough the sensor 10. As outlined above, sucrose molecules, whenpresent, will be trapped in the porous silicon layer 12. This results ina change in the optical properties of the layer 12 manifest as a changein the value of n_(avg)d, in this instance an increase. When sucrose isremoved from its environment, and pure water is again pumped through thesensor 10, the trapped sucrose molecules are effectively flushed out andthe value of n_(avg)d returns to its original value. Strictly, the valueof n_(avg)d falls at time t₂ to the projected average background level74, which in fact is lower than at t₁. This reduction in background is along term drift in the observed level of n_(avg)d it occurs because theporous silicon structure is slightly unstable and is either slowlyetched away or oxidised when in solution.

FIG. 5 is a plot similar to that shown in FIG. 4 but taken from anactual set of results obtained using the chemical sensor andwater/sucrose solution as described above and analysed using the priorart FFT technique. It gives graphs of optical thickness n_(avg)d (shownas nd) against time: with the equipment used the interferogram scan ratewas linear in wavenumber, so the time axis is also a measure ofwavenumber. Similar remarks apply to other drawings with time axes.

The results 80 take the form of three separate traces 82, 84, 86corresponding to pumping sucrose solutions through the sensor film 12,the solutions being 0.1% (upper trace 82), 0.05% (middle trace 84) and0.03% (lower trace 86) sucrose in water. Points t₃ and t₄ are marked onthe time axis as indicators of the times at which the sucrose solutionswere first introduced and then removed from the sensor environmentrespectively. FIG. 5 shows that it is not possible to detect reliablyfrom the traces 82, 84, 86 when either of the environmentalchanges—sucrose introduction or removal—occurred.

As noted previously, there are in fact several defects in theapplication of the FFT analysis method to data obtained from thebiosensor 10 of FIG. 1. Data is collected at regular wavelengthintervals, which inevitably corresponds to a non-linear series ofwavenumbers Analysis therefore requires interpolation of data betweenregular sampling points in the Fourier transformed spectrum.Interpolation assumes a particular correlation function (e.g. a linearrelationship between adjacent samples), which is imposed on the data andwhich will differ from the true correlation function. Such filteringinevitably distorts the shape of the transformed spectrum. In order toachieve the desired sensitivity in measurement of n_(avg)d to 1 part in10 000, the position of the desired peak in the FFT spectrum must bedetermined with an accuracy better than 1/10000th of the width of aFourier resolution cell. This is very difficult to achieve byinterpolating between the points in the FFT spectrum. In addition, theobserved interference patterns are not exactly sinusoidal, primarily dueto the fact that the refractive index varies with wavelength but alsodue to some systematic noise in the measurements. The FFT method hasbeen found to be sensitive to this distortion and also to noise. It isclear from FIG. 5 that this prior art technique fails to extractsufficiently clean results from measured interferometric spectra toallow transduction of the molecular trapping, in this instance, to beobserved.

Finally, although a simple monolayer biosensor system with only atwo-beam reflected interference pattern is described herein, there is noobvious extension of the prior art FFT method to the analysis of morecomplicated interference patterns obtained from more complex structures,e.g. a multilayer porous silicon structures. These structures give riseto interference patterns which can be far from sinusoidal andaccordingly use of a FFT would simply not be of any assistance in theiranalysis.

By way of contrast to the prior art, the present invention does not relyon interpolation of a Fourier transform spectrum. Nor even necessarilyon the Fourier transform technique at all, although in approximatelyperiodic systems there are advantages in making use of it to someextent. There are two stages involved in implementing this invention,and novel aspects of both stages will be described.

In the mathematical analysis which follows, the expression ‘function’means a closed form continuous or analytic function, and a set of values(perhaps measurements, results or intermediate results) is representedas a vector such as X having members (vector elements) X_(j). Such avector may be used to replace the function, and represent itdiscontinuously as a set of discrete points such as X_(j). The term‘entity’ is used to mean an item which is a vector member, a vector or afunction.

In general, the method of analysis of the invention assumes that a setof discrete observations or measurements is made of an observed variablex, the set being represented by the data vector X. The observed variablex may be regarded as a member of the set X or equivalently as an elementof the vector X, and it is derived by measurements upon a physicalsystem characterised by an unknown and variable parameter ω. Each memberof X represents a measurement made at a respective value of aninterrogation or measurement parameter of the kind k, so that the jthmember of X, X_(j)=x(ω, k_(j)). In accordance with the invention, it isassumed that the function x, or x(ω, k) since it is a function of ω andk, can be expressed as the sum of an imperfectly known function ƒ(ω, k)and a random noise process η(k), i.e.:x(ω, k)=ƒ(ω, k)+η(k)  (2)

Without affecting the generality of the method of the invention, it isapparent that the detector 24 in the apparatus of FIG. 1 detects arespective voltage at each wavenumber. In this example the parameterX_(j) may therefore be designated as a detector voltage measured at aparticular value of k, which is designated as the wavenumber of thatmeasurement. Upper case X is a vector having vector elements or memberswhich collectively form a set of voltages providing individual points onand used herein to define a whole interferogram.

The interferogram 26 therefore constitutes a set X of individualvoltages X_(j), and each voltage X_(j) is measured at a respectivediscrete wavenumber in a wavenumber set. The parameter ω may thereforebe designated as the optical thickness of the porous silicon layer 12i.e. n_(avg)d when its pores are filled with a solvent. Noise is alsoincluded in the method as the parameter η(k), which is independent of ω.

It is an important characteristic of the invention that it makes avirtue out of necessity. The invention is directed to detectingperturbations to a physical system manifesting themselves as very smallchanges to some physical parameter, optical thickness in the presentexample. The invention exploits the fact that the parameter changes tobe detected are very small by using that as a basic assumption inanalysing experimental data. In this connection, for small changes in ωabout a reference value ω₀, ƒ(ω,k) can be expanded about ω₀ and Equation(2) can be rewritten:x(ω,k)=ƒ(ω₀ ,k)+(ω−ω₀)ƒ′(ω₀ ,k)+η(k)+O((ω−ω₀)²)  (3)where ƒ′(ω₀,k) is the derivative of ƒ(ω,k) with respect to ω at ω=ω₀ andO((ω−ω₀)²) denotes terms of the second and higher orders in (ω−ω₀), i.e.(ω−ω₀)², (ω−ω₀)³ etc. For small (ω−ω₀), the terms O((ω−ω₀)²) can beignored. If both ƒ(ω₀,k) and ƒ′(ω₀,k) can be estimated, then an estimateof the parameter change or perturbation (ω−ω₀) between ω and ω₀ isobtainable from Equation (3) by evaluating:

$\begin{matrix}{\left( {\omega - \omega_{0}} \right)_{e} = {\frac{\left( {{x\left( {\omega,k} \right)} - {f\left( {\omega_{0},k} \right)}_{e}} \right)}{{f^{\prime}\left( {\omega_{0},k} \right)}_{e}}\mspace{14mu}{for}\mspace{14mu}{any}\mspace{14mu}{k.}}} & (4)\end{matrix}$where the subscript e denotes an estimated value.

In order to implement the model represented by Equation (4), two stagesare employed, stages 1 and 2. Stage 1 is concerned with finding areference value ω₀ with respect to which the parameter of interest ωwill be determined as a difference (ω−ω₀). Stage 2 relates to evaluatingEquation (4) above or an equivalent at specific values of the parameterk: this is in order to determine the parameter change or shift (ω−ω₀)from ω₀ caused by a change in this example to the optical properties ofthe porous layer 12. Initially a vector F is used to denote a set ofdiscrete values or points on the function ƒ(ω₀,k)_(e) at values of kdefined by k=k_(j), j=1 to N_(k). Similarly, the vector F′ denotes theset of values of the function ƒ′(ω₀,k)_(e), at the same set of values ofk. To implement stages 1 and 2, two series of data vectors X, i.e. twoseparate series of interferograms, are obtained: the first such seriesis a series of calibration data vectors, X^(j) _(c), j=1 to N_(c), whichdefines a background state of the sensor system. The second such seriesis a series of test data vectors, X^(j) _(t), j=1 to N_(t), in which aperturbation manifest as a change in a system parameter is to bemeasured. Here N_(c) is the number of calibration data vectorscollected, and N_(t) is the number of test data vectors collected.Typical values in the examples described herein are N_(c)=20, andN_(t)=60. Each interferogram has 1050 points, so each data vector X^(j)_(t) has 1050 members X^(j) _(tl).

Various alternatives are available to provide a reference value ω₀: twosuch will be described herein with reference to the chemical sensorsystem of FIG. 1. Clearly many alternative approaches may be used, beingsomewhat dependent on the nature of the system under analysis and theinformation available with which to formulate an estimate.

FIG. 6 is a flow chart of an example of steps involved in making theestimates F and F′ in applying the present invention to the poroussilicon biosensor system of FIG. 1. A first approach to stage 1 of theanalysis is represented by steps 90 to 96 of this drawing. Firstly, at90 a set of calibration data (CD) vectors X^(j) _(c), j=1 to N_(c), iscollected using the apparatus of FIG. 1. Each vector X^(j) _(c) ismeasured at a different time, and the noise in each measurement (seeEquation (2)) will be different. These vectors are obtained by scanningthrough the measurement wavenumber range several times and recordingdetector voltages X_(j) at specific wavenumbers. They are collectedwhile pure water flows through the biosensor 10, and thus are known tobe due to a hydrated porous silicon layer 12. Test data vectors X^(j)_(t), j=1 to N_(t) taken in biosensor environments other than pure watermay be collected at this time. In fact they are ideally collected afterone set of calibration data vectors and before another, to provide atwo-sided calibration with intervening test. It is found however, that aone-sided calibration (i.e. that performed using data taken eitherbefore or after the test measurements) is also effective. In eithercase, the test data vectors are not subjected to analysis until stage 2of this analysis system.

In general application of the invention, it is apparent that the actualnature of the calibration data vectors will depend on the system underinvestigation.

The next step 92 is calculation of an average calibration data vector(CD_(av)) from the collected set of CD vectors or interferograms of thekind X^(j) _(c). This is carried out by adding the voltages X₁ at afirst measurement wavenumber k₁ on all the calibration interferogramsand dividing by the total number. This is then repeated for all otherwave numbers. An average interference pattern CD_(av) or averagecalibration data vector results from this step, and provides a basis forderiving an estimate vector F of the function ƒ(ω,k) at ω=ω₀. F has arespective member at each value of k_(j), j=1 to N_(k).

In analysing the signal from the biosensor system detector 24, thenature of the relationship between variables makes it convenient to usea Fourier transform function in order to derive an estimate of thereference value ω₀. Thus a fast Fourier transform (FFT) spectrum of theaverage CD_(av) vector is computed at 94 (other spectral analysistechniques may also be used). In this average FFT spectrum, the positionand phase φ of its largest amplitude peak is determined at 96. Thisposition is used to calculate the reference value ω₀, which is theoptical thickness of the porous layer 12 when pervaded by a puresolvent, in this case water.

A second approach to stage 1 comprises analysing biosensor signals in anovel manner. It obviates the need for interpolation of an FFT: this isbecause in certain circumstances it has been found that it isunnecessary to Fourier transform an entire data spectrum. A Fouriertransform can instead be calculated directly in a region of interest.This novel analysis component will be referred to herein as a directFourier transform.

In application to the biosensor 10, this approach recognises that theapproximate optical thickness of the porous silicon layer 12 will beknown beforehand from the layer thickness, porosity and solvent type.From this known thickness an expected Fourier transform peak positioncan be deduced. The Fourier transform peak is not expected to shiftsignificantly in response either to an analyte becoming bound to arecognition agent or a solute molecule becoming trapped in the poroussilicon layer 12, because the resulting change in optical propertieswill normally be too small for this. Accordingly, the Fourier transformcan be computed directly for a series of optical thickness values in apredetermined range of such values: the range is centred on the opticalthickness value which corresponds to the expected location of theFourier transform peak. Estimates of the peak position are fine-tuned byuse of a line search or iterative peak finding position. An importantadvantage of this example is that it does not use interpolation. Insteadthe Fourier transform is computed at a set of data points in the regionof interest directly from the collected data, i.e. at data pointsinstead of between them so interpolation error does not arise.

Each interference pattern detected by the biosensor system of FIG. 1 andused in the derivation of the traces shown in FIG. 5 has been analysedusing the above-described direct Fourier transform method. That is, foreach data vector, the peak position (and hence value of n_(avg)d) hasbeen extracted using this method of analysis of the invention instead ofthe prior art FFT technique. FIGS. 7 and 8 display the results of thisanalysis: FIG. 7 illustrates three traces 87 a, 88 a and 89 a ofn_(avg)d, again corresponding to immersion of the sensor 10 in purewater (background calibration over times 1 to 16 and 45 to 60) and invarious sucrose solutions (measurements over times ˜20 to ˜45). For theplots 87 a, 88 a and 89 a, the sucrose solutions are 0.1%, 0.05% and0.03% sucrose in water respectively. FIG. 8 illustrates three traces 87b, 88 b, 89 b obtained by dividing each of 87 a to 89 a of FIG. 7 by arespective mean background value, i.e. in each case the mean of thevalues obtained over times 1 to 16 and 45 to 60. It can be seen fromFIGS. 7 and 8 that a distinct improvement is immediately apparent overthe prior art analysis system used to generate the results shown in FIG.5. Moreover, it is further clear that in the cases of all three solutionstrengths, a threshold can be set which will distinguish between waterand sucrose environments. It is to be noted that this improvement waseffected by implementing stage 1 of this invention only.

Note that it is not essential to use a Fourier transform at this stage.All that is required is a method of determining the initial referencevalue ω₀. Thus signals from other perturbed physical systems may beamenable to other forms of analysis. Two embodiments of stage 1 havebeen described herein however in relation to an approximately periodicinterference signal, and accordingly Fourier transforms have been usedfor convenience.

The second, novel, stage of this analysis technique arises fromrecognition that in many systems, it is not necessary to measure anabsolute value of a parameter for every configuration: instead it isenough to establish an approximate reference or background value andthen to measure variations about this reference. That is, this stage isapplicable whenever variations in an observed parameter can be seen asperturbations relative to a reference value of that parameter. It isaccordingly to be noted that it is relevant to many perturbationdetection applications, and not merely to analysis of roughly periodicsignals. For convenience however, roughly periodic signals will continueto be used as illustrative examples herein.

Using this analysis technique, parameter variations are extracteddirectly. It is therefore unnecessary to calculate the Fouriertransform, provided that some other way is available to generate areference signal. In particular, it is not necessary to make theassumption that the signal being detected is a pure sine wave, thatassumption being inherent in the Fourier transform technique. This leadsto an improved signal to noise ratio in the processed signal andconsequently, better detectability of event thresholds.

There will now be described the steps involved in deriving the estimatesF and F′ for use in deriving an error vector E with a jth member E_(j).This employs Equation (5), which is obtained by evaluating Equation (4)at each value of k=k_(j), j=1 to N_(k):

$\begin{matrix}{E_{j} = {\left( {\omega - \omega_{0}} \right)_{e} = {\frac{\left( {X_{tj} - F_{j}} \right)}{F_{j}^{\prime}}.}}} & (5)\end{matrix}$

Ideally the estimates F and F′ (with jth values F_(j) and F_(j)′) areobtained from a number of data vectors (interference pattern valuesets): they are then weighted according to their signal to noise ratiosand averaged, although it is believed that such weighting would make thesystem more sensitive to the uncertainty in knowing ƒ′(ω₀,k)_(e).

Referring once more to FIG. 6, and as described earlier, by Step 96 avalue and phase φ for the reference value ω₀ has been derived in stage 1of the analysis technique. The range of optical thicknesses expected tobe encountered using the biosensor 10 is known, and therefore the rangeof interferogram peak positions expected is also known. A range ofdiscrete values of ω spanning this thickness range and centred on ω₀ isselected at 98 and a corresponding set S of sinusoidal data vectors isgenerated. For this purpose parameters of the kind δω₁ are defined whichare offsets expressing the range of discrete values of ω as individualoffsets from ω₀. This provides a set of sinusoidal functionss(ω_(l),k_(j)) given by:{s(ω_(l) ,k)}={sin((ω₀+δω_(l))k+φ)}  (6)where {} indicates a set, δω_(l) is the lth in a set of offsets{δω_(l)}, each offset may be positive or negative and represents arelatively (< 1/10) small departure from ω₀, and φ is the phase of thereference value, as referred to above. The corresponding set S ofsinusoidal data vectors is obtained by evaluation of these functionss(ω_(l),k_(j)) at each value of k=k_(j), j=1 to N_(k) according toEquation (7)S={S ^(l) }:S ^(l) _(j) ={s(ω_(l) ,k _(j))}  (7)

Once more, the selection of sine waves is specific to the nature of theinterference pattern for which there is sought a reference functionestimate, ƒ(ω₀,k)_(e) evaluated at each value of k_(j), j=1 to N_(k) togive a reference vector F. The biosensor 10 is assumed to give rise toan interference pattern which is an approximate cosine function, so areference vector for such a pattern will have the form of a set ofdiscrete points on such a function. The method of the invention detectsperturbations with respect to the reference vector, i.e. smalldepartures from the reference vector: such departures will be associatedwith the reference derivative vector, sine functions in this example.

At step 100 the calibration data is filtered using the set S. That is,each observation or measurement vector X^(j) _(c), j=1 to N_(c), made incollecting the calibration data vectors (interference patterns), isassumed to be a respective linear combination of all members of the setS. This treats each X^(j) _(c), j=1 to N_(c), as being projected on toeach sinusoidal data vector in the set S acting as a basis vector ineach case. Components in the direction of each sinusoidal data vector inS are therefore computed for each observation vector X^(j) _(c), j=1 toN_(c). This step 100 produces filtered calibration data (CD^(f)) whichagain is averaged as described earlier at a following step 102 toproduce an estimate of the imperfectly known reference data vector, F interms of the reference value ω₀ and the set of offsets {δω_(l)} Thefiltering removes explicit dependence upon the interrogation parameteror wavenumber k: that is, the observation vector X has been transformedfrom k space to a projected observation vector Z in the space of thesinusoidal data vectors, with their dependence on selected (discrete)values of δω, according to Equation (8):

$\begin{matrix}{Z_{l} = {{z\left( {\omega,{\delta\;\omega_{l}}} \right)} = {\sum\limits_{j}{{x\left( {\omega,k_{j}} \right)}{\sin\left( {{\left( {\omega_{0} + {\delta\;\omega_{l}}} \right)k_{j}} + \phi} \right)}}}}} & (8)\end{matrix}$

Equation (8) represents the lth member Z_(l) of the projectedobservation vector Z as taking the value z(ω,δω_(l)) at δω=δω_(l); eachvalue z(ω,δω_(l)) (general value l) is in turn a respective summation

$\sum\limits_{j}{{x\left( {\omega,k_{j}} \right)}{\sin\left( {{\left( {\omega_{0} + {\delta\;\omega_{l}}} \right)k_{j}} + \phi} \right)}}$over products of a (jth) member x(ω,k_(j))=X_(j) of the observationvector X with a respective (lth) member sin((ω₀+δω_(l))k_(j)+φ) of theset S of sinusoidal data vectors. Each projected observation vector Zcorresponds to a respective observation vector X.

At step 102 the projected observation vectors Z corresponding in themanner described above to the calibration data (CD) constitute thefiltered calibration data (CD^(f)). Then F is computed by adding uprespective members of each projected observation vector Z^(j) _(c), anddividing the total by the number of such vectors, N_(c), according toEquation (9):

$\begin{matrix}{F_{l} = {\sum\limits_{j = 1}^{N_{c}}{{z^{j}\left( {\omega,{\delta\;\omega_{l}}} \right)}/N_{c}}}} & (9)\end{matrix}$where again F_(l) denotes the lth member of the data vector F, andz^(j)(ω,δω_(l)) is the lth member of the jth projected calibration datavector Z^(j) _(c), evaluated at δω=δω_(l). Thus Equation (9) establishesF for use as the reference vector in Equation (5). The reference vectorF, a series of discrete values derived for each value of at δω=δω_(j),j=1 to 2m+1, is now used to represent the estimated reference functionƒ(ω₀,δω)_(e), a continuous function.

FIG. 9 is a flow chart of process steps for direct extraction of thevalue of the parameter ω from successive test data vectors X^(j) _(t),j=1 to N_(t), once the reference vector F has been established followinga route similar to that outlined with reference to FIG. 6. In FIG. 9, ata first step 110, test or ‘real’ data RD is acquired in the same way asdescribed earlier for calibration data CD, except that real data isobtained with a small amount of sucrose present in the otherwise purewater solvent. The real date RD is filtered using the set of sinusoidaldata vectors S, as defined in Equation (7). That is, each measurement isagain assumed to comprise a linear combination of members of the set S.The set S does not contain high and low frequency terms, and so bothhigh and low frequency noise are suppressed in the filtered real data.Subsequent processing is performed on real data RD^(f) filtered in thisway. The filtering step 110 is not essential but it improves signal tonoise ratio. At step 112, the reference vector F is subtracted from eachdata vector Z^(j) _(t), for each for each j, j=1 to N_(t) of thefiltered real data RD^(f): this subtraction provides a difference vectorwhich is the error vector E^(j), for each j, j=1 to N_(t), i,e. arespective error vector is derived for each test interferogram. FromEquation (5), it can be seen that each member of the error vector isproportional to (ω−ω₀).

It is now assumed that (before projection on to {s(ω_(l),k)}, and neverexplicitly computed) the estimated reference function ƒ(ω₀,k)_(e) is ofthe form cos(ωk). That this is reasonable approximation is illustratedin FIG. 10, where the estimated reference function ƒ(ω₀,k)_(e) isplotted as a function of wavenumber k (solid line) compared to a cos(ωk)approximation (dashed line). The derivative of cos(ωk) with respect toω, i.e.

${\frac{\partial}{\partial\omega}{\cos\left( {\omega\; k} \right)}},$is k sin(ω₀k) when evaluated at the reference point ω₀. This estimatedderivative function is denoted ƒ′(ω₀,k)_(e), and when evaluated at eachvalue of k=k_(j), j=1 to N_(k) yields the derivative vector F′.

In order to make a proper interpretation of small differences between aprojected observation vector Z and a projected reference vector F, thederivative vector F′ is also projected at step 114 onto the set ofsinusoidal data vectors S of Equation (7) to give members of F′ asfollows:

$\begin{matrix}{F_{l}^{\prime} = {\sum\limits_{j}{{f_{j}^{\prime}\left( {\omega_{0},k_{j}} \right)}_{e}{\sin\left( {{\left( {\omega_{0} + {\delta\;\omega_{l}}} \right)k_{j}} + \phi} \right)}}}} & (10)\end{matrix}$where again F_(l)′ denotes the lth member of the derivative vector F′.

Referring now to FIG. 11, the vector F′ represents the projectedestimated derivative function ƒ′(ω₀,δω)_(e), evaluated at each value ofδω=δω_(j), j=1 to 2m+1. The function ƒ′(ω₀,δω)_(e) represented by F′ isshown as a solid curve 140 plotted as a function of ω₀+δω: it is asinc(ω₀+δω)-like function, i.e. it is of the form (sin θ)/θ. It iscompared to similar plots 142 (dots), 144 (dashes) and 146 (dash-dot)for three error functions of the form e(ω,δω), which when evaluated ateach value of δω=δω_(j), j=1 to 2m+1, are represented by the errorvectors E. The plots 142, 144 and 146 are taken respectively from apreceding calibration interval (interval from times 1 to 16 in FIG. 8)the measurement interval (times 20 to 40 in FIG. 8) and a secondcalibration interval (times 48 to 60). The plot 140 has a peak 140 ashown in FIG. 12, in which like located regions 142 a to 146 a are alsoshown. The error functions e(ω,δω) are equal to the numerator ofEquation (4), i.e. to (x(ω,k)−ƒ(ω₀,k)_(e)), and from that equation theyare also equal to ƒ′(ω₀,k)_(e)(ω−ω₀)_(e). From this it can also be seenthat if the functions ƒ′(ω₀,k)_(e), x(ω₀,k)_(e) and ƒ(ω₀,k)_(e) are allprojected on to the sine function space {s(ω_(l),k)} removing thedependency on wavenumber k, then:e(ω,δω)=ƒ′(ω₀,δω)_(e)(ω−ω₀)_(e) =x(ω,δω)−ƒ(ω ₀,δω)_(e)  (11)is true for all values of δω. This means one can also combine resultsfor all values of δω by integrating Equation (11). Thus, integrating anddividing both sides of Equation (11) by ∫ƒ′(ω₀,δω)_(e):

$\begin{matrix}{\left( {\omega - \omega_{0}} \right)_{e} = {\frac{\int{e\left( {\omega,{\delta\;\omega}} \right)}}{\int{f^{\prime}\left( {\omega_{0},{\delta\;\omega}} \right)}_{e}} = \frac{\int\left( {{x\left( {\omega,{\delta\;\omega}} \right)} - {f\left( {\omega_{0},{\delta\;\omega}} \right)}_{e}} \right.}{\int{f^{\prime}\left( {\omega_{0},{\delta\;\omega}} \right)}_{e}}}} & (12)\end{matrix}$

Moreover, summation over a subset of the δω values in theoriginally-selected range of interest (that surrounding ω₀) will alsohold true. It is not essential to carry out a summation, any member ofan error vector can be used to provide an estimate (ω−ω₀)_(e), butsumming over a number of selected low noise members improves accuracyfor the usual statistical reasons. Thus, converting continuous functionsin Equation (12) to data vectors and integrals to summations over theirmembers:

$\begin{matrix}{\left( {\omega - \omega_{0}} \right)_{e} = {\frac{\sum E}{\sum F^{\prime}} = \frac{{\sum X} - F}{\sum F^{\prime}}}} & \left( {12a} \right)\end{matrix}$where the summations are of vector members corresponding to a subset ofδω values chosen as follows: the chosen subset comprises those δω valuesfor which the respective associated values of the error function arelarge relative to residual noise. Error function values are large overthe peak 140 a of the sinc function ƒ′(ω₀,δω)_(e) at 140. Therefore, atstep 116 (FIG. 9) a summation is made of values of the error functione(ω,δω) at values of δω that span the sinc peak 140 a. This correspondsto summing appropriate members of the error vector E for aninterferogram, since, as previously explained, each member of E, E_(j)is equal to the value of the function e(ω,δω) evaluated at the value ofδω=δω_(j), j=1 to 2m+1.

In this example of the invention, an interferogram and the data vectorXwhich represents it has 1050 measurement points/members, and errorvector members corresponding to 101 these points are selected forsummation. The sinc peak 140 a occurs at about 4650 nm corresponding toω₀. It is about 700 nm wide, and requiring 101 points spread overapproximately the upper half of the peak 350 nm wide (full width halfmaximum) gives δω values with a spacing of 3.5 nm. The general value ofδω, i.e. δω_(j), is then given to be a wavelength increment or decrementof (3.5(j−1)−175) nm. The maximum value of δω is approximately ω₀/27; δωis therefore no more than relatively small compared to ω₀, less than onetenth of it.

The summation over the sinc peak 140 a provides the numerator of theright hand side of Equation (10). It gives an indication that ω haschanged, but not the absolute value of that change. For many purposessuch an indication is all that is needed, such as for example when analarm signal is required to indicate a change from a required state hasoccurred.

It is an important property of the invention that signal to noise ratiocan be improved as described by only using relatively large values oferror functions and discarding those with poorer signal to noise ratios.

To provide an absolute value of the shift in ω using Equation (12a),rather than just an indication that it has changed, the normalisationfactor ΣF′ in the denominator of the right hand side of that equation isrequired. At step 118, ΣF′ is evaluated and applied as follows. It isequivalent to evaluating the function ƒ′(ω₀,δω)_(e) for δω in the chosensubset of δω values referred to above, these function values beingsummed. This corresponds to summing appropriate members of F′, since, aswas previously explained, each member of F′, F_(j)′ is equal to thevalue of ƒ′(ω₀,δω)_(e) evaluated at δω=δω_(j), j=1 to 2m+1. Thesummation is performed over the subset of δω values over which membersE_(j) of the error vector E were summed. The sum of the error vectorvalues (numerator on the right hand side of Equation (12a)) is thendivided by the resulting normalisation factor ΣF′ to give (ω−ω₀), andhence ω itself as ω₀ is known.

There are a number of enhancements which may be added to the analysisdescribed above. The first enhancement concerns random amplitudefluctuations which appear in the interference patterns. See, forexample, the pattern shown in

FIG. 2. Referring to FIGS. 6 and 9 once more, in this enhancement thesubtraction at step 112 is replaced by a projection. That is, instead ofderiving the error function e(ω,δω) by subtracting ƒ(ω₀,δω)_(e) fromeach filtered interference pattern using Equation (11), each filteredpattern is projected onto a space orthogonal to that of F, the referencevector. This calculation is equivalent to normalising each projecteddata vector Z by dividing by an estimate of the unknown amplitude ofeach projected data vector Z and then subtracting F.

Secondly, in place of the static reference function ƒ(ω₀,δω)_(e), atsteps 102 and 112 a linearly time-varying reference functionƒ(ω₀,δω,t)_(e) is used, where t is time. This will be described in termsof functions for convenience, although strictly speaking it should be interms of data vectors as described earlier for the time-invariant case.It is achieved by assuming the calibration data follows a best-fit slopeover time, and this is taken into account at the averaging step 92. Thetime-varying reference function or estimated interference patternƒ(ω₀,δω,t)_(e) is consequently defined to be the sum of the staticreference function ƒ(ω₀,δω)_(e) computed from Equation (9), and alinearly time-varying pattern v(ω₀,δω)_(e) calculated from Equation(13).

$\begin{matrix}{{v\left( {\omega_{0},{\delta\;\omega}} \right)}_{e} = \frac{\left( {\sum\limits_{j = 1}^{N_{c}}{\left( {j - {\left( {1 + N_{c}} \right)/2}} \right)\left( {{z^{j}\left( {\omega,{\delta\;\omega}} \right)} - {f\left( {\omega_{0},{\delta\;\omega}} \right)}_{e}} \right)}} \right)}{\left( {\sum\limits_{j = 1}^{N_{c}}\left( {j - {\left( {1 + N_{c}} \right)/2}} \right)^{2}} \right)}} & (13)\end{matrix}$

In Equation (13), N_(c) is the number of observation data vectors X_(c)obtained in calibration. This Equation shows that the functionv(ω₀,δω)_(e) is obtained by subtracting ƒ(ω₀,δω)_(e) from z^(j)(ω,δω),which is the member of the jth projected calibration data vector Z^(j)corresponding to δω, for each value of j from 1 to N_(c), scaling theeach resulting difference by multiplying by the current value ofj−(1+N_(c))/2, and then adding up all the scaled differences. Theresulting total is divided by the sum of the squares of the values ofj−(1+N_(c))/2. Equation (14) is then used to compute ƒ(ω₀,δω,t)_(e):ƒ(ω₀ ,δω,t)_(e)=ƒ(ω₀,δω)_(e) +tv(ω₀,δω)_(e)  (14)

Computer software implementing this embodiment of the invention (i.e.including both the enhancements described above) was developed using“Matlab®” produced by Mathworks Inc., an American corporation. Thissoftware is attached at Annex A. It clearly demonstrates that equationsgiven in the foregoing description can be evaluated by an appropriatecomputer program recorded on a carrier medium and running on aconventional computer system. Such a program is straightforward for askilled programmer to implement without requiring invention, because thecode exemplifies the equations and the equations themselves are wellknown computational procedures.

FIG. 13 is an illustration of spectra similar to those plotted in FIGS.5, 7 and 8, but now obtained using the above-described enhanced analysistechnique. This drawing shows three plots 160, 162 and 164 of the kinddescribed earlier: i.e. optical thickness is plotted against time(wavenumber) for three different sucrose solutions with compositionspreviously given. The plots 160, 162 and 164 exhibit far more pronounceddiscontinuities when the sensor is exposed to or washed free of therespective sucrose solutions (see times 15-20 and 42-45) than theequivalent in FIGS. 5, 6 and 8.

Although described herein in relation to detection of events using abiosensor, it is to be understood that perturbation detection accordingto this invention is susceptible to far wider application. It shouldprove useful in any application in which it is necessary to detect smallchanges in the response of a sensor and where conventionalhigh-resolution methods fail because of uncertain calibration. Onepossible application is to laser vibrometry.

It will be apparent to one skilled in the art that the algorithm used inthis analysis system can be modified for use with interference (orother) patterns which are not simple sinusoids. All that is needed is toobtain an estimate of a state of a system in terms of a parameter ofinterest and its derivative with respect to that parameter. Thisprovides ƒ(ω₀,k)_(e) and ƒ′(ω₀,k)_(e) or equivalents thereof and setsthe projection space, and then the same strategy applies.

Referring now to FIG. 14, there are shown two interference patternsindicated generally by 200 (solid curve) and 202 (chain curve)respectively. The patterns 200 and 202 show reflectivity plotted againstwavelength in nm. The pattern 200 was measured as described withreference to FIG. 1, except that the single porous silicon layer 12 wasreplaced by a sensor film of porous silicon having eighteen layers andconfigured as a Bragg mirror. Each layer consisted of first and secondfilms with different porosities, so the sensor film was effectivelythirty-six films or layers in which porosity alternated between twovalues in successive layers. In other respects the sensor film wasilluminated with light and the interference pattern 200 was detected asdescribed earlier. The pattern 200 was obtained while pure water waspassed through the sensor film: it provides a reference spectrum and wasderived by averaging data in calibration mode.

The pattern 202 was calculated as a reflectivity predicted by a simpleanalytic model for the sensor film referred to above. For thismultilayer sensor film the sensed parameter is a wavelength λ₀ at whicha peak of a spectrum occurs, and an interrogation parameter is providedby illuminating light wavelength λ incident on the film. Although itdoes not provide a perfect fit to pattern 200, the model pattern 202 issufficiently accurate for use to provide an acceptable approximation toƒ′(λ₀,λ)_(e), which serves the same purpose as ƒ′(ω₀,k)_(e) in Equation(4). The model pattern 202 is ƒ(λ₀,λ)_(e), and it can be shown that itis described by Equation (15) below:ƒ(λ₀,λ)=K ² sin h ²(y)/(Γ² cos h ²(y)+δ² sin h ²(y))  (15)where K=2(n1−n2)/λ,

-   -   δ=(π/L)(2(d₁n₁−d₂n₂)/λ−1)    -   Γ=√ K²−δ²)    -   L=d₁+d₂, is the thickness of each layer (film thicknesses d₁ and        d₂ having different refractive indices n₁ and n₂),    -   d₁=λ₀/(4n₁), is the thickness of the first film (quarter        wavelength at λ₀),    -   d₂=λ₀/(4n₂), is the thickness of the second film (quarter        wavelength at λ₀),    -   n₁ is the refractive index of the first film,    -   n₂ is the refractive index of the second film,    -   y=mLΓ, and    -   m is the number of layers in the multilayer sensor, i.e. m=18 in        this example.

By differentiating the model pattern 202 with respect to λ at λ₀,ƒ′(λ₀,λ)_(e) is obtained as defined in Equation (16)

$\begin{matrix}{{f^{\prime}\left( {\lambda_{0},\lambda} \right)} = \frac{\left( {4{\pi\delta}\; K^{2}{{\sinh(y)}/\lambda}} \right)\left( {{\sinh(y)} - {{mL}\;{{\Gamma cosh}(y)}}} \right)}{\left( {{\Gamma^{2}{\cosh^{2}(y)}} + {\delta^{2}{\sinh^{2}(y)}}} \right)^{2}}} & (16)\end{matrix}$

Referring to FIG. 6 once more, the multilayer sensor data is processedin the same way as that from a single layer sensor, except that insteadof steps 94 and 96 the wavelength λ₀ corresponding to the peak of theaverage interference pattern is computed. In step 98 a set P of basisvectors is computed by evaluating a set of functions of the formp_(l)(λ₀−δλ_(l),λ)=ƒ′(λ₀−δλ_(l),λ) _(e) at discrete values of λ andδλ_(l), where δλ_(l) is the lth in a series of small wavelength changes.In step 100 the calibration data is filtered using P, in the mannerdescribed in Equation (10), where p_(l)(λ₀−δλ_(l),λ_(j)) replacessin((ω₀+δω_(l))k_(j)+φ). In step 102 ƒ(λ₀,δλ)_(e) is computed byprojection instead of ƒ(ω₀,δω)_(e).

Likewise in FIG. 9 at step 110 the basis set P is used instead of S. Atstep 112 ƒ(λ₀,δλ)_(e) is subtracted instead of ƒ(ω₀,δω)_(e), to givee(λ₀,δλ). At step 114 ƒ′(λ₀,λ)_(e) is projected on to P to obtainƒ′(λ₀,δλ)_(e), which is sinc-like. In step 116 components of e(λ₀,δλ)are summed, and in step 118 the result is normalised to provide anestimate of λ₀+δλ, indicating the required shift.

In this example a model pattern was calculated as for the pattern 202 toestimate responses to levels of sucrose dissolved in water.

FIG. 15 is a graph of response (arbitrary units) of the multilayersensor processing system described above against time sample number(equivalent to time and wavenumber). Peaks 170 and 172 indicate theresponse of this system to different concentrations of sucrose in water,0.1% and 0.05% respectively. The use of a multilayer porous silicon filmis advantageous in this application as the reflectivity and consequentlythe signal to noise ratio is improved as compared to the single layerfilm, allowing a more sensitive measurement system to be constructed.

Annex A MATLAB code % Script to read in psi experiment data; generate atime-varying estimate % of the mean water spectrum and subtract it fromeach snapshot vector % in turn. Then correct for amplitude fluctuations.% NB graphics removed % read in data from files % 0.1% sucrose%filename=‘data010.txt’ %outfile=‘output010.txt’ %n=60;%plot_title=‘0.1% sucrose’; % 0.05% sucrose %filename=‘data005.txt’%outfile=‘output005.txt’ %n=58; %plot_title=‘0.05% sucrose’; % 0.03%sucrose filename=‘data003.txt’ outfile=‘output003.txt’ n=60;plot_title=‘0.03% sucrose’; % assume the data has two regions ofcalibration data w1=17; % last calibration sample in first block w2=50;% first calibration sample in second block w3=n; % last calibrationsample in second block % computation options op1=0; % one sidedcalibration; op1=0 => two sided calibration op2=1; % computetime-varying reference; op2=0 => fixed referenceXa=read_data(filename,n+1); % separate the data D=Xa(:,2:n+1); % fromthe wavelength (nm) info, and convert to wavenumber in nm{circumflexover ( )}−1 t=1.0./Xa(:,1); t=t−min(t); % since the waveforms we areinterested in are cosine waves % find the frequency which fits best fromthe specified range % lower frequency (nm) fl=2000.0; % upper frequency(nm) fu=15000.0; % step size (nm) df=10.0; freq=[fl:df:fu]; % find bestfit frequency and starting phase C=exp(j*2*pi.*t*freq); d=mean(D.‘).’; %average spectra vs=C‘*(d−mean(d)); % remove mean [mr,kk]=max(abs(vs)); %find frequency with maximum amplitude phase=angle(vs(kk)); % estimatephase % set upper & lower limits for integrationfln=freq(kk)−1/(2*max(t)); fun=freq(kk)+1/(2*max(t)); % matrix Sprojects onto space of sine waves S=imag(C*exp(j*phase)); Dp=S‘*D; %obtain linear time varying average if op1 == 0 g0=[[1:w1],[w2:w3]]; elseg0=[1:w1]; end offset=mean(g0); g=g0−offset; % X is the calibration dataonly X=Dp(:,g0); vm=mean(X.‘).’; [m,n]=size(X); Y=zeros(size(X)); %project out mean for k=1:n Y(:,k)=X(:,k)*vm.‘*vm/(vm.’*X(:,k))−vm; end %calculate linear time-varying component Z=repmat(g,[m,1]).*Y;vd=mean(Z.‘).’/mean(g.{circumflex over ( )}2); M=zeros(size(Dp)); % nowproject off a time-varying reference for each data vector fork=1:size(Dp,2) v=vm+op2*(k-offset)*vd;M(:,k)=Dp(:,k)*v.‘*v/(v.’*Dp(:,k))−v; end % reduce range of interest forintegration q=[fix((fln−fl)/df)+1:fix((fun−fl)/df)+1];s=S‘*(S(:,kk).*t); es=−length(t)*sum(M(q,:))/(2*pi*df*mr*sum(s(q))); %plot output figure plot((es).{circumflex over ( )}l) xlabel(‘sampleno.’) ylabel(‘response’) title(plot_title) % write to filefid=fopen(outfile,‘w’); fprintf(fid,‘%6.3f\n’,es) fclose(fid)

1. A method of detecting perturbation of a physical system from areference state associated with a reference parameter (ω₀) to aperturbed state associated with a perturbed parameter (ω), the methodincluding the steps of: a) deriving the reference parameter (ω₀); b)deriving a reference entity (F) which describes the system's state atthe reference parameter (ω₀); c) deriving an error entity (E) from thedifference between the reference entity (F) and a measurement-relatedentity (Z) associated with a perturbed state of the system; and d) usingthe error entity (E) to provide an indication of whether or not thesystem state has been perturbed from the state associated with thereference parameter (ω₀), with the perturbed parameter (ω) having becomeunequal to the reference parameter (ω₀).
 2. A method according to claim1 wherein the error entity (E) is a vector with multiple members (E_(j),j=1 to N_(k)) and step d) includes calculating the sum of the errorvector's members (E_(j)) to obtain an indication of whether or not theperturbed parameter (ω) has become unequal to the reference parameter(ω₀).
 3. A method according to claim 2 wherein the error entity (E) hasmembers characterised by relatively high signal to noise ratio comparedother possible members which might otherwise be selected to derive it.4. A method according to claim 3 wherein the error entity (E) hasmembers derived from a region of a peak in a derivative (F′) of thereference entity (F).
 5. A method according to claim 1 wherein prior toderivation of the error entity (E) in step c), the measurement-relatedentity (Z) is normalised by projection on to a space orthogonal to thatof the reference entity (F).
 6. A method according to claim 1 whereinstep d) includes determining the difference (ω−ω₀) between the perturbedand reference parameters (ω, ω₀) by error entity normalisation withrespect to an entity which is a summation of elements of an entity (F′)representing a derivative (ƒ′(ω₀,δω)_(e)) of a reference function(ƒ(ω₀,δω)_(e)) evaluated at the reference parameter (ω₀), the referencefunction (ƒ(ω₀,δω)_(e)) being represented by the reference entity (F).7. A method according to claim 1 wherein step b) incorporates; a)predicting from the reference parameter (ω₀) a position of a peak in aFourier transform spectrum of observation data; b) selecting a subset ofthe observation data over the predicted peak and calculating a directFourier transform (as herein defined) for the subset; and c) analysingthe direct Fourier transform in order to derive a more accuratedetermination of the position of the peak.
 8. A method according toclaim 1 wherein the reference entity (F) is derived in step b) by aprocess which includes filtering by projection of an entity on to a setS of pre-determined entities having a range of arguments all of whichdiffer from the reference parameter (ω₀) by less than one tenth of thereference parameter (ω₀).
 9. A method according to claim 8 wherein stepc) of claim 1 comprising calculating an error entity (E) includesderiving the measurement-related entity (Z) by a process which includesdata filtering by projection of a measurement entity (X) on to the set Sof pre-determined entities.
 10. A method according to claim 1 whereinthe reference entity (F) is derived in step b) from an average of aseries of observation data vectors (X).
 11. A method according to claim10 wherein the observation data vectors (X) are time-varying and this istaken into account when their average is derived.
 12. A method accordingto claim 11 wherein the reference entity (F) is a sum of a staticreference entity and a linearly time-varying entity.
 13. A methodaccording to claim 1 wherein the system is a sensor (10) pervaded by amedium having variable composition.
 14. A method according to claim 13wherein the sensor is a porous silicon sensor (10) with pores pervadedby the medium which is a solvent, the perturbed and reference parameters(ω, ω₀) are optical thicknesses of a region of the sensor (10) pervadedby different solvent compositions and in which interference patterns aregenerated for derivation of the reference entity (F) and error entity(E).
 15. Apparatus for detecting perturbation of a physical system froma reference state associated with a reference parameter (ω₀) to aperturbed state associated with a perturbed parameter (ω), the apparatusincluding: a) a monitoring device for monitoring the physical system toobtain indications of the system's state; b) a computer systemprogrammed to process the system state indications to: i) derive thereference parameter (ω₀); ii) derive a reference entity (F) whichdescribes the system's state at the reference parameter (ω₀); iii)derive an error entity (E) from the difference between the referenceentity (F) and a measurement-related entity (Z) associated with aperturbed state of the system; and iv) use the error entity (E) toprovide an indication of whether or not the system state has beenperturbed from the state associated with the reference parameter (ω₀),with the perturbed parameter (ω) having become unequal to the referenceparameter (ω₀).
 16. Apparatus according to claim 15 wherein the errorentity (E) is a vector with multiple members (E_(j), j=1 to N_(k)) andthe computer system is programmed to calculate the sum of the errorvector's members (E_(j)) to obtain an indication of whether or not theperturbed parameter (ω) has become unequal to the reference parameter(ω).
 17. Apparatus according to claim 16 wherein the error entity (E)has members characterised by relatively high signal to noise ratiocompared other possible members which might otherwise be selected toderive it.
 18. Apparatus according to claim 17 wherein the error entity(E) has members derived from a region of a peak in a derivative (F′) ofthe reference entity (F).
 19. Apparatus according to claim 15 whereinthe computer system is programmed to normalise the measurement-relatedentity (Z) by projection on to a space orthogonal to that of thereference entity (F) prior to derivation of the error entity (F). 20.Apparatus according to claim 15 wherein the computer system isprogrammed to determine the difference (ω−ω₀) between the perturbed andreference parameters (ω, ω₀) by error entity normalisation with respectto an entity which is a summation of elements of an entity (F′)representing a derivative (ƒ′(ω₀,δω)_(e)) of a reference function(ƒ(ω₀,δω)_(e)) evaluated at the reference parameter (ω₀), the referencefunction (ƒ(ω₀,δω)_(e)) being represented by the reference entity (F).21. Apparatus according to claim 15 characterised in the computer systemis programmed to derive the reference entity (F) by; a) predicting fromthe reference parameter (ω₀) a position of a peak in a Fourier transformspectrum of observation data; b) selecting a subset of the observationdata over the predicted peak and calculating a direct Fourier transform(as herein defined) for the subset; and c) analysing the direct Fouriertransform in order to derive a more accurate determination of theposition of the peak.
 22. Apparatus according to claim 15 wherein thecomputer system is programmed to derive the reference entity (F) by aprocess which includes filtering by projection of an entity on to a setS of pre-determined entities having a range of arguments all of whichdiffer from the reference parameter (ω₀) by less than one tenth of thereference parameter (ω₀).
 23. Apparatus according to claim 22characterised in the computer system is programmed to calculate theerror entity (E) by deriving the measurement-related entity (Z) by aprocess which includes data filtering by projection of a measuremententity (X) on to the set S of pre-determined entities.
 24. Apparatusaccording to claim 15 wherein the computer system is programmed toderive the reference entity (F) from an average of a series ofobservation data vectors (X).
 25. Apparatus according to claim 24wherein the observation data vectors (X) are time-varying and thecomputer system is programmed to take this into account when derivingtheir average.
 26. Apparatus according to claim 25 wherein the referenceentity (F) is a sum of a static reference entity and a linearlytime-varying entity.
 27. Apparatus according to claim 15 wherein thesystem is a sensor (10) pervaded by a medium having variablecomposition.
 28. Apparatus according to claim 27 wherein the sensor is aporous silicon sensor (10) with pores pervaded by the medium which is asolvent, the perturbed and reference parameters (ω, ω₀) are opticalthicknesses of a region of the sensor (10) pervaded by different solventcompositions and in which interference patterns are generated forderivation of the reference entity (F) and error entity (E).
 29. Acomputer program product comprising a computer readable mediumcontaining computer readable instructions for use in detectingperturbation of a physical system from a reference state associated witha reference parameter (ω₀) to a perturbed state associated with aperturbed parameter (ω), the computer readable instructions being forcontrolling computer apparatus to implement the steps of: a) derivingthe reference parameter (ω₀); b) deriving a reference entity (F) whichdescribes the system's state at the reference parameter (ω₀); c)deriving an error entity (E) from the difference between the referenceentity (F) and a measurement-related entity (Z) associated with aperturbed state of the system; and d) using the error entity (E) toprovide an indication of whether or not the system state has beenperturbed from the state associated with the reference parameter (ω₀),with the perturbed parameter (ω) having become unequal to the referenceparameter (ω₀).
 30. A computer program product according to claim 29wherein the error entity (E) is a vector with multiple members (E_(j),j=1 to N_(k)) and step d) includes calculating the sum of the errorvector's members (E_(j)) to obtain an indication of whether or not theperturbed parameter (ω) has become unequal to the reference parameter(ω₀).
 31. A computer program product according to claim 30 wherein theerror entity (E) has members characterised by relatively high signal tonoise ratio compared other possible members which might otherwise beselected for deriving it.
 32. A computer program product according toclaim 31 wherein the error entity (E) has members derived from a regionof a peak in a derivative (F′) of the reference entity (F).
 33. Acomputer program product according to claim 29 wherein prior toderivation of the error entity (E) in step c), the measurement-relatedentity (Z) is normalised by projection on to a space orthogonal to thatof the reference entity (F).
 34. A computer program product according toclaim 29 wherein step d) includes determining the difference (ω−ω₀)between the perturbed and reference parameters (ω, ω₀) by error entitynormalisation with respect to an entity which is a summation of elementsof an entity (F′) representing a derivative (ƒ′(ω₀,δω)_(e)) of areference function (ƒ(107 ₀,δω)_(e)) evaluated at the referenceparameter (ω₀), the reference function (ƒ(ω₀,δω)_(e)) being representedby the reference entity (F).
 35. A computer program product according toclaim 29 wherein step b) incorporates: a) predicting from the referenceparameter (ω₀) a position of a peak in a Fourier transform spectrum ofobservation data; b) selecting a subset of the observation data over thepredicted peak and calculating a direct Fourier transform (as hereindefined) for the subset; and c) analysing the direct Fourier transformin order to derive a more accurate determination of the position of thepeak.
 36. A computer program product according to claim 29 wherein stepb)'s deriving of the reference entity (F) includes filtering byprojection of an entity on to a set S of pre-determined entities havinga range of arguments all of which differ from the reference parameter(ω₀) by less than one tenth of the reference parameter (ω₀).
 37. Acomputer program product according to claim 36 wherein the step ofcalculating the error entity (E) includes deriving themeasurement-related entity (Z) by a process which includes datafiltering by projection of a measurement entity (X) on to the set S ofpre-determined entities.
 38. A computer program product according toclaim 29 wherein the step of deriving the reference entity (F) employsan average of a series of observation data vectors (X).
 39. A computerprogram product according to claim 38 wherein the observation datavectors (X) are time-varying and this is taken into account when theiraverage is derived.
 40. A computer program product according to claim 39wherein the reference entity (F) is a sum of a static reference entityand a linearly time-varying entity.
 41. A computer program productaccording to claim 29 wherein the system is a sensor (10) pervaded by amedium having variable composition.
 42. A computer program productaccording to claim 41 wherein the sensor is a porous silicon sensor (10)with pores pervaded by the medium which is a solvent, the perturbed andreference parameters (ω, ω₀) are optical thicknesses of a region of thesensor (10) pervaded by different solvent compositions and in whichinterference patterns are generated for derivation of the referenceentity (F) and error entity (E).